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Abstract. 

Atmospheric flows are governed by the equations of fluid dynamics. These 
equations are nonlinear, and consequently the hierarchy of cumulant equations is 
not closed. But because atmospheric flows are inhomogeneous and anisotropic, the 
nonlinearity may manifest itself only weakly through interactions of non-trivial 
mean fields with disturbances such as thermals or eddies. In such situations, 
truncations of the hierarchy of cumulant equations hold promise as a closure 
strategy. 

Here we show how truncations at second order can be used to model and 
elucidate the dynamics of turbulent atmospheric flows. Two examples are 
considered. First, we study the growth of a dry convective boundary layer, which 
is heated from below, leading to turbulent upward energy transport and growth of 
the boundary layer. We demonstrate that a quasilinear truncation of the equations 
of motion, in which interactions of disturbances among each other are neglected 
but interactions with mean fields are taken into account, can capture the growth of 
the convective boundary layer. However, it does not capture important turbulent 
transport terms in the turbulence kinetic energy budget. Second, we study the 
evolution of two-dimensional large-scale waves, which are representative of waves 
seen in Earth’s upper atmosphere. We demonstrate that a cumulant expansion 
truncated at second order (CE2) can capture the evolution of such waves and 
their nonlinear interaction with the mean fiow in some circumstances, for example, 
when the wave amplitude is small enough or the planetary rotation rate is large 
enough. However, CE2 fails to capture the flow evolution when strongly nonlinear 
eddy-eddy interactions that generate small-scale filaments in surf zones around 
critical layers become important. Higher-order closures can capture these missing 
interactions. 

The results point to new ways in which the dynamics of turbulent boundary 
layers may be represented in climate models, and they illustrate different classes 
of nonlinear processes that can control wave dissipation and angular momentum 
fluxes in the upper troposphere. 


Keywords: turbulence, closure, quasi-linear approximation, atmospheric boundary 
layer, atmospheric convection, large-scale atmospheric circulation, jets 
Submitted to: New J. Phys. 



2 


1. Introduction 


Atmospheric flows shape Earth’s climate and are governed by the equations of 
fluid dynamics, the Navier-Stokes equations augmented by the Coriolis force and 
thermodynamic equations [e.g., Ooyama 2001 Vallis 2006 Pauluis 2008 , and 


equations for the microphysical processes describing, for example, the formation and 
re-evaporation of cloud droplets Pruppacher et al., 1998 . They span an enormous 


range of length scales, from the micrometers of droplet formation to the planetary 
scale. Temporal variations range from microseconds at the smallest scales to tens of 


years on the largest scales [e.g., Klein 2010 . Atmospheric processes are tightly coupled 
across all of these scales. For example, cloud droplets scatter sunlight and absorb 
infrared radiation, thereby affecting Earth’s radiative balance globally; conversely, 
planetary-scale dynamics affect where and how clouds form. Current climate models 
cannot resolve all relevant scales. They resort to the direct simulation of dynamics 
on scales of tens of kilometers and larger, while representing smaller-scale processes 
such as turbulence in clouds and boundary layers semi-empirically [e.g., Beljaars 


[TM^[U arratt[ |1994[ |Smith[ |1997[ [Lappen and Randall[ |2001[ |Soares[ |2004[ |Siebesrna] 


et al., 2007 . However, the larger-scale dynamics of weather systems, with timescales 


of minutes, are simulated explicitly even when only their longer-term statistics—the 
climate—is ultimately of interest. 

Two scientific objectives would be beneficial to achieve. First, it would be 
desirable to obtain more accurate and more physically motivated models of the 
interactions between the larger scales that can currently be resolved in climate models 
and the smaller scales that cannot be resolved. Inaccuracies in how these interactions, 
in particular in clouds and boundary layers, are represented in climate models are the 


largest source of uncertainties in climate projections [e.g., Stephens 2005 Bony et al. 


|2006[ [Soden and Held[ |2006[ [Webb et al.[ |2013[ [Stevens and Bon;^ |2013[ [Vial et aTj 
2013 Brient et al. 2016 . Improving the representation of such interactions would 


have an enormous societal benefit. Second, it would be desirable to devise ways of 
calculating climate statistics more directly, rather than through the direct simulation 
of weather systems and accumulation of their statistics, as is currently done. This may 
in the long run lead to faster climate simulations. In the shorter term, it may lead to 
insight into how climate is maintained and how it varies on timescales of seasons to 
millennia. 

Both objectives require the development of closure models for the hierarchy 
of statistical moment or cumulant equations associated with the equations of fluid 
dynamics. This hierarchy is in principle infinite because of the quadratic nonlinearity 
of the Navier-Stokes equations. Numerous ways of closing the hierarchy of moment 
or cumulant equations in a variety of circumstances have been proposed [see, e.g., 
Frisch 1995 Pope 2000 Lesieur 2008 . Many of them concern flows that are 


assumed statistically homogeneous and isotropic, as an idealized benchmark from 
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which the development of closures for more realistic applications can proceed [e.g., 
[Orszag 1970, 1973 . However, closures for homogeneous and isotropic turbulence 


often may not be easier to obtain than closures for more realistic flows: Mean fields 
in homogeneous and isotropic turbulence can, without loss of generality, be taken to 
be zero; only higher-order statistics of the flows are of interest. Isotropic turbulence 
cannot equilibrate with any imposed driving and dissipation through interaction with 
mean flows; rather, it must equilibrate through nonlinear interactions across scales. 
By contrast, turbulence in the atmosphere usually interacts strongly with non-trivial 
mean fields, which include, for example, atmospheric jet streams or the thermal 
stratification of the atmosphere. Because interactions between turbulent fluctuations 
and non-trivial mean fields have the potential to be important, many atmospheric 
flows may be less strongly nonlinear than the oft-studied prototype problems of three- 


dimensional turbulence [e.g.. 

Pedlosky 

1970 

Farrell 

1987 

Farrell and loannou 

1993 

Randel and Held| 

1991 

Schneider and Walker 

2006 . Moreover, already the mean 


fields (e.g., mean temperatures and winds) are of primary interest for understanding 
climate, though, of course, higher moments (e.g., temperature extremes) also remain 
important to understand. 

Because the nonlinearity of turbulent interactions in many atmospheric flows 
may be limited, truncating the hierarchy of moment or cumulant equations at a low 
order has potential to be successful. Here we explore the feasibility of truncations at 
second order—that is, neglecting third-order nonlinearities in second-order covariance 
equations—in two prototype problems of atmospheric flows. The first is a turbulent 
convective boundary layer, with scales of motion on the order of meters to a kilometer. 
The second is a model of large-scale turbulence in the upper atmosphere, with 
scales of motion on the order of hundreds to thousands of kilometers. These two 
problems involve disparate phenomenologies and force balances. For example, the 
boundary layer can be taken to be unaffected by the planetary rotation, whereas the 
planetary rotation and Coriolis forces are fundamental for the large-scale turbulence 
in the upper atmosphere. Yet the problems share that turbulent fluctuations interact 
strongly with a non-trivial mean state—a thermal stratification in the first case, and 
an atmospheric jet in the second case. Because of the strength of this interaction, 
truncations of moment equations at second order already achieve some success in 
capturing the statistics of these flows. It is essential in these truncations that nonlocal 
and anisotropic covariation of turbulent quantities (e.g., in waves or convective plumes) 
are retained. Such nonlocal truncation of the moment or cumulant equations at 


second order is known as second-order cumulant expansion (CE2) [Marston et al. 

2008 

Marston 

2012 

Srinivasan and Young 

2012 

Tobias and Marston 2013 or 

stochastic structural stability Theory (S3T) [Farrell and loannou 

2003 

2007 

Bakas 

and loannou ^ 

>014 

Constantinou et al.[ 2014a . CE2 and S3T differ in that S3T 


attempts to represent missing eddy-eddy interactions, whereas CE2 sets them to zero. 
CE2 is a realizable closure in that its equations are the exact moment equations 
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of a realizable system, the quasi-linear (QL) system that corresponds to the original 
equations of motion. The QL approximation retains the interaction of turbulent 
fluctuations with a mean flow but neglects the interactions of turbulent fluctuations 
among each others. CE2 and S3T were successful in explaining some aspects of zonal 


jet dynamics in rotating flows [e.g.. 

Farrell and loannou 

2003 

to 

o 

o 

Srinivasan and 

Young 

2012 

Tobias and Marston 

2013 

Bakas and loannou 

2013 , without relying 

on eddy-eddy interactions and inverse cascades Rhines 

1975 

Vallis and Maltrud 


1993 . QL approximations of large-scale atmospheric dynamics, sometimes with added 


damping and stochastic forcing, were partially successful in reproducing aspects of the 


atmospheric climate and its variability [e.g., 

Whitaker and Sardeshmukh 1998 

Zhang 

and Held 

1999 

DelSole 

2001 

O’Gorman and Schneider 

2007 . At smaller scales, QL 


approximations also capture sheared stably stratified flows when the dynamics involve 


the linear excitation and absorption of internal gravity waves Orr , 1907 Lindzen 1988 


Bakas and loannou, 2007 . They also reproduce aspects of thermal convection, such as 


the dependence of the heat flux on the Rayleigh number [e.g., Malkus 1954 Herring 


1963 Toomre et al. 1977 Busse 1978 Niemela et al. 2000 


In what follows, we derive the CE2 closure for Boussinesq flows, present fully 
nonlinear and QL simulations of a dry convective boundary layer using large-eddy 
simulations (LES), and study the evolution of a large-scale wave disturbance on a 
zonal jet representative of the upper troposphere. The results will demonstrate the 
potential and limitations of CE2 approaches. 


2. Cumulant expansion of Boussinesq flow 

2.1. Boussinesq flow 

Atmospheric flows have low Mach number, so sound waves are generally unimportant 
for the dynamics. It is therefore common to study atmospheric dynamics with 
approximations to the Navier-Stokes equations that filter sound waves. The simplest 
such approximation, which ignores all density variations except where they affect the 
buoyancy of air masses, is the Boussinesq approximation [Boussinesq [1872] . The 
Boussinesq equations are obtained by expressing density p{r,t) = pQ + Sp{r,t) as 
a sum of a constant density po in a reference state and fluctuations 6p(r, f) about 
it, assuming pressure variations in the reference state are hydrostatically balanced, 
and retaining only the leading-order terms in density and pressure fluctuations in an 


expansion of the Navier-Stokes equations. The resulting equations are [e.g., Vallis 


2006 


clu 


-I- (u • V)u + 2Q, X u = —V<I> -I- 5k -)- 

(momentum equation) 

(la) 

V • u = 0 

(continuity equation) 

(lb) 

db 

— +VL-Vb = Fb 
at 

(thermodynamic equation) 

(Ic) 
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Here, u denotes the three-dimensional velocity, Sp the pressure perturbation associated 
with the density perturbation Sp, $ = Sp/po the potential of the pressure-gradient 
accelerations, and b = —gSp/po the buoyancy {g is an effective gravitational 
acceleration and k is the local vertical). The reference frame rotates with the constant 
angular velocity Cl of the planetary rotation, as a result of which Coriolis accelerations 
(217 X u) appear in the momentum equation. Centrifugal accelerations are subsumed in 
g, the effective gravitational acceleration. The terms and on the right-hand side 
represent dissipation and forcing terms (for example, friction and diabatic heating). 
The continuity equation reduces to the condition that the flow u is non-divergent. 
Sound waves are filtered from these equations because no time derivative appears in 
the continuity equation. In effect, the speed of sound is taken to be infinite, so that 
pressure adjusts instantaneously across the flow domain: it can be determined from a 
Poisson equation obtained by taking the divergence of the momentum equation and 
using the non-divergence condition to eliminate the time derivative. 

To write the equations in a synthetic way, we introduce the state vector 

’3/ = (u , u , w , b)^ (2) 

of the flow that contains all prognostic variables in Cartesian coordinates (the 
superscript (•)^ indicates the transpose). The set of equations Q can then be written 
compactly as 

5^-fV-(4'(g)u)=L'®'-V$-fF (3a) 

at 

V • u = 0, (3b) 


where the outer product (o) of the two vectors \I/ and u is defined as 


(g) u = 


• • • 5'iM3 

'I'4Ui • • • 'I'4M3 


(d>iUj)i<i<4. 

l<i<3 


( 4 ) 


The components of the divergence of the second-order term V • (tl' (g) u) are 

[V • (’®' 0 u)], = A , (5) 

where summation over repeated indices is implied. This notation extends the usual 
advection operator of a scalar field by a non-divergent flow to the advection of a 
vector held. The linear operator L contains accelerations owing to the Coriolis force, 
buoyancy, as well as nonconservative terms (e.g., friction or diabatic heating) that are 
linear in the state vector. The vector F contains all other nonconservative terms. We 
have expanded u, the outer product and the V • in Cartesian coordinates for simplicity, 
however the formalism is coordinate independent. 
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Despite their relative simplicity, the Boussinesq equations are commonly used 
to study turbulence in boundary layers, in which density variations are weak. They 
also underlie classical conceptual models of large-scale atmospheric dynamics (e.g., 
quasigeostrophic models), in which they become quantitatively inaccurate but still 
qualitatively capture important atmospheric phenomena such as Rossby waves and 
baroclinic instability. The Boussinesq equations are well suited for our exposition 
of cumulant expansion approaches because they capture the essential nonlinearity of 
atmospheric flows: the conservative quadratic nonlinearity of the advection operators 
(u • V)u and u • V6. 


2.2. Averaging operator 

Our interest is not in individual details of the atmospheric flows under consideration 
but in their statistics, including mean values and higher moments. Therefore, we 
define an averaging operator, denoted with an overline (•), and decompose any scalar 
field /(r, t) into a mean and a fluctuating part, 

f{r,t) = f{r,t)+ (6) 


The mean is in general a function of space and time. The fluctuating part is commonly 
referred to as an eddy. The averaging operator is taken to satisfy, for all scalar fields 


/(r, t) and g{r,t) and any constant c, the Reynolds properties Monin and Yaglom 

[IM] : 


c = c (7a) 

cf + g = cf + g (linearity) (7b) 

fg = fg (7c) 

9f = df (commutation with derivatives) (7d) 


Properties c) imply that the averaging operator is a projection operator and so is 
idempotent / = /. They make it possible to define the Reynolds decomposition 


fg = fg + fg'- 


( 8 ) 


For a vector quantity ’®', the average ’5' is the component-wise average. 

The choice of average is unspecified as long as it satisfies the Reynolds properties. 
Conceptually, ensemble averages are statistically meaningful, and they naturally 
satisfy the Reynolds properties. In practice, however, they can be difficult to obtain. 
Averages over sufficiently long times in statistically stationary (or slowly varying) flows 
or over sufficiently large regions in statistically homogeneous flows are more commonly 
used in practice, and also approximately satisfy the Reynolds properties. In concrete 
calculations, we choose the averages that are natural given the statistical symmetries of 
the problem under consideration. For example, in flows that are statistically invariant 
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under translations along a spatial coordinate (e.g., along latitude circles), an average 
along that spatial coordinate suggests itself. 

For more general flow equations with a variable density, the averaging operator 
above has to be replaced by a density-weighted average, to obtain consistent equations 
of motion for the statistical moments that resemble the Boussinesq moment equations 
formally. An example for the anelastic approximation is provided in [Appendix A[ 

2.3. Cumulant expansion 

First cumulant The first cumulant is the mean ’®'(r,t), for which the equations of 
motion are obtained by averaging the equations of motion (§: 

- _ - - — 

-f V • (’*' (g) u) = -V • (g) u') L ^ - V0 F, (9a) 

V • u = 0. (9b) 

This involves the covariance V• (’®'' (g) u'), which arises from the quadratic nonlinearity 
of the equations of motion. It represents eddy fluxes, for example, arising from 
advection of momentum fluctuations by the eddies (fluctuations) themselves. Because 
the equation for the mean involves a covariance, it is not closed. 


Second cumulant The second cumulant is the second central moment, or the 
covariance 

C(ri, r 2 , t) = ’®"(ri, t) (g ^'(r 2 , t). (10) 

We only consider the prognostic variables here, because the diagnostic variables (and 
hence their statistics) can be obtained from them. We also only consider equal-time 
cumulants, that is, equal-time covariances between prognostic variables at the two 
points ri and r 2 , which need not be equal. The covariance tensor is symmetric, 

C(ri,r2,t) = C^(r2,ri,t). (11) 

We additionally define the auxiliary covariances, 

C*(ri,r 2 ,t) = $'(ri,t)^'(r 2 ,t) (12a) 

C“(ri,r 2 ,t) = ’®"(ri,f) g u'(r 2 ,t). (12b) 


The first, C^, contains additional information about covariation of the prognostic 
fields '5' with the pressure potential $. These covariances can be calculated from C 
and 'h because the pressure potential is a diagnostic variable in the Boussinesq 
approximation. The second, C“, represents covariation of the velocity field with other 
prognostic variables and is already contained in C. 

The second cumulant equation can be obtained from the equations of motion ^ 
by evaluating 




C(ri,r2,t) 


'i"(ri,t) 


94''(r2,t) 


94''(ri, t) 


^'(r2,t). 


dt 


dt 


(13) 













Discarding terms that are third order in fluctuating quantities, one obtains 


— C(ri,r 2 ,t) + Vri • [C(ri,r 2 ,t) 0 u(ri,t)] = -C“(ri,r 2 ,t) [V§(r 2 ,t)]'^ 

+ LC(ri,r 2 ,t) + VriC'^(ri,r 2 ,t) + F'(ri,t) (g) '4"(r2,t) + {ri <—^ r 2 }, (14) 

where {ri <—>■ r 2 } indicates the terms obtained by interchanging ri and r 2 , which 
are necessary to ensure the symmetry 0 of the covariance tensor. The third-order 
term C(ri,r 2 ,t) 0 u(ri,f) is defined as in Cartesian coordinates, 

[C(ri,r2,f) g u(ri,t)]p ,j . = Cp,,(ri,r2, t)Mi(ri, t), (15) 

and its divergence by 

{Vri • [C(ri, r 2 , t) g u(ri, f)]}p_^ = ^ [Cp,qUi\ . (16) 

The flux C(ri, r 2 , t) g u(ri, t) represents the transport of spatial eddy correlations by 
the mean flow at ri. The term — C“(ri,r 2 ,t) • V [^(r 2 ,t)]^ represents generation of 
covariance C(ri,r 2 ,t) by advection down mean-flow gradients. 

Additionally, continuity &) implies that 

Vr, • C“(ri,r 2 , t) = Vr, • C“(r 2 ,ri,t) = 0. (17) 

The set of equations for the first and second cumulants in this form is closed because 
the third cumulants, which would ordinarily appear in the second cumulant equations 
owing to the quadratic nonlinearity of the equations of motion, have been discarded. 
This second-order truncation of the otherwise infinite hierarchy of cumulant equations 
is referred to as a second-order cumulant expansion (CE2). In that CE2 assumes the 
first and second cumulants suffice to describe the flow statistics, it makes a normal 
approximation to the hierarchy of equations for the flow statistics. 

CE2 equations To summarize, the CE2 equations, with the second cumulant 
substituted in ([^), are given by 

+ V • ['S'(r, t) g u(r, t)] = -V • C"(r, r, t) 

-fL'5'(r,t)-V$(r,t)-hr(r,<) (18a) 

V • u = 0 (18b) 

^C(ri,r 2 ,<) -I-Vri • [C(ri,r 2 ,t) gu(ri,t)] = -C“(ri,r 2 ,t) [V’l'(r 2 , t)] 

-fLC(ri,r2,t) -b VriC'^(ri,r2,t)-l-F'(ri,t) g ’S"(r2,t) -t- {ri <—^ r2} 

(18c) 

Vr, • C“(ri,r2,t) = Vr, • C“(r2,ri,t) = 0. (18d) 
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This set of equations involves the 16 terms in C and the eight covariances C* (ri, r 2 , t) 
and C^(r 2 ,ri,t). The 16 components of the second cumulant C are prognostic 


(evoloving according eq. (18c))), and the other 8 covariances components involve 


diagnostic correlations between the pressure potential and each of the other fields. The 
8 corresponding diagnostic Poisson equations are obtained by taking the divergence 
with respect to ri or r 2 of the equation for C“(ri,r 2 ,t) that can be extracted from 


(18c), and in which time tendencies vanish because of the non-divergence constraint 


(18d). 


Physical properties CE2 is a realizable approximation, in the sense that the implied 
probability distribution functions are positive Marston et al. 2014 . The first 


cumulant equations (18 1 , b) are unchanged from the fully nonlinear system. Therefore, 


first-order invariants (mass and momentum) are conserved by the CE2 system in the 
absence of nonconservative effects. The third cumulants, which are neglected in the 


second cumulant equations (18:, d), redistribute second-order invariants (e.g., energy) 
among scales but do not generate or dissipate these invariants. Therefore, second-order 
invariants such as energy are likewise conserved by the CE2 system in the absence of 


nonconservative effects [see Marston et al. 2014 . 


2.4- Quasi-linear approximation 

The CE2 equations can also be obtained by directly approximating the original 
equations of motion (§, making what has come to be called the quasi-linear 


(QL) approximation [O’Gorman and Schneider 

2007 

Srinivasan and Young 

2012 

[Constantinou et al. 2014b 

Ait-Chaalal and Schneider 

2015 . The QL approximation 


keeps nonlinear terms in the equation ([^ for the mean 'i'. However it linearizes the 
equation for the eddies obtained by substracting ([^) from ([^), 


_ _ 

—-bV-(’®'(g)u')+V-(4''(g)u)-h[V-(^'0u')-V-(’4"(g)uOl = L(19) 
ot 

by setting 

V • (’4''(g) u') = V • (’S"(g)u'). (20) 

Hence, the QL approximation amounts to replacing in the Reynolds decomposition of 
the nonlinear term, 

4'(g)u = ’4'(g)u-|-^(g)u'-|-’4''(g)u-|-’4''(g)u', (21) 

the eddy-eddy interaction (g) u' by its mean effect \I/' ® u'. Under the QL 
approximation, the Boussinesq equations (§ can then be written as 


_ 

— -b V • (’4' (g) u) -H V • (’4" (g u' - '4'' (g u') = L ’4' - -b F, (22a) 
ot 

V • u = 0. (22b) 
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Because the QL equations retain as the only nonlinear interaction the interaction 
between eddies and mean fields, the corresponding cumulant equations are closed at 
second order, meaning no third-order terms appear in the second-order equations. 


The first two cumulant equations are exactly the CE2 equations (18). This makes it 


possible to simulate flows whose statistics satisfy the CE2 equations (181 simply by 


simulating the QL equations (22). 

The QL truncation does not necessarily imply that eddy-eddy interactions and 
third-order correlations are equal to zero. However their evolution is decoupled 
from that of lower-order cumulants. This has to be kept in mind when interpreting 
instantaneous fields and statistics of flows simulated by the QL equations. 


2.5. Differences between CE2 and S3T 

Stochastic Structural Stability Theory (S3T) is a second-order statistical approach to 
turbulent flows that is closely related to CE2. CE2 and S3T are sometimes presented as 
being equivalent in the literature because they share a similar mathematical formalism. 

However, CE2 and S3T differ in that S3T includes a small-scale stochastic forcing 
that is white in time and represents the scattering by missing eddy-eddy interactions. 
The resulting additional energy injection is balanced by large-scale linear damping. 
The stochastic forcing allows one to define rigorously an ensemble average over its 
realizations and permits a semi-analytical treatment of the second-order equations, 
whose solutions depend on the existence and statistics of the stochastic forcing. 

By contrast, CE2 uses the same forcing in the truncated equations as in the fully 
nonlinear equations, without attempting to parameterize eddy-eddy interactions. This 
choice is made for two reasons. First, a stochastic noise is not necessarily a realistic 
model for eddy-eddy interactions because these interactions do not inject energy but 
only redistribute it among scales. Second, forcing the flow at small scales might 
appear unnatural for a wide range of planetary flows, which are forced on larger 
scales. For example, large-scale motion in Earth’s atmosphere is essentially driven 
by the planetary-scale radiative imbalance between the equator and the poles, rather 
than by energy injection at small scales. 

Another difference is that S3T is generally applied to flows for which the linear 
operator representing eddy-mean flow interactions does not have unstable modes (the 
mathematical development of S3T relies on non-normal growth and decay of stable 
modes, excited by the stochastic noise). In this framework, S3T has been used in 


barotropic rotating flows to study instabilities of zonal flows Bakas and loannou 2013 


Parker and Krommes 2014 and the dynamics of zonal [Farrell and loannou 2003 


2007 Constantinou et al. 2014a and non-zonal Bakas and loannou 2014 coherent 


structures. Nevertheless, the simulation of unstable flows at second order is possible 
and well-defined mathematically. Hence, CE2 (or its QL analogue) has been applied 
to unstable flows in barotropic Marston et al. 2008 and baroclinic jO’ Gorman and 







































Schneider 2007 Ait-Chaalal and Schneider 2015 settings. Evaluating the ability of 


CE2 in capturing unstable flows is at the cornerstone of the present paper. 


3. Dry convective boundary layer 

The dynamics of boundary layers and clouds involve flow scales of order meters to 
a few kilometers. In addition, condensate formation and evaporation take place at 
microscales. By contrast, typical climate models have a horizontal resolution of order 
100 km. Therefore, the dynamics of boundary layers and clouds are subgrid-scale 
processes in climate models, which must be represented parametrically in terms of 
the resolved large-scale dynamics [e.g., Beljaars 1992 Smith, 1997 . Uncertainties 


about these parameterization schemes are the dominant contributor to uncertainties 


in climate change projections [e.g.. 

Stephens 

2005 

Bony et al. 

2006 

So den and Held 

2006 

Webb et al. 

2013 

Stevens and Bony| 2013 '' 

/ial et al. 

2013| Brient et al.| 

Mel. 


Most current parameterization schemes for the dynamics of clouds and boundary 
layers truncate the hierarchy of moment (or cumulant) equations at first order 
and represent the second-order subgrid-scale fluxes appearing in the first-order 
equations semi-empirically. For example, turbulent fluxes in boundary layers are often 
represented as diffusive fluxes down mean gradients, with diffusivities estimated from 
approximate spatially local second-order equations for turbulence kinetic energy [e.g., 
Mellor and Yamada 1982 Beljaars 1992 . Turbulent fluxes in convective clouds, on 


the other hand, are often represented as vertically non-local entraining plumes that 
extend deep into the boundary layer. The parameters determining the mass fluxes in 


the plumes are estimated from energy equations [e.g., Arakawa and Schubert 1974 


Gregory 1997 


CE2 may offer a path toward improved and unified representations of subgrid- 
scale turbulent fluxes (e.g. in boundary layers and in convective clouds) in climate 
models, because CE2 explicitly retains spatial nonlocality and interactions between 
fluctuations (plumes or eddies) and the environment (mean field). Here we compare 
a QL simulation and a fully nonlinear simulation of the simplest convective boundary 
layer, a dry convective boundary layer, and demonstrate the potential and limitations 
of CE2 to represent the statistics of its dynamics. We conducted a large-eddy 
simulation (LES) of a dry convective boundary layer growing into a stable background 


stratification as a heat flux is imposed at the surface Soares 2004 . We then compared 


this fully nonlinear simulation with a simulation in which the equations were replaced 


by the corresponding QL equations (22). 


3.1. Large-eddy simulations 

Setup The LES code solves the anelastic equations and is described in|Pressel et al.| 


2015 . Like the Boussinesq approximation, the anelastic equations are non-hydrostatic 


and filter sound waves by neglecting dynamic density variations except where they 
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affect buoyancy. In contrast to the Boussinesq approximation, the background state 
depends on the vertical coordinate. Hence, the anelastic equations can capture the 
dynamics of flows with substantial stratification and so are better suited to study 
atmospheric convection. The anelastic approximation and its cumulant expansion are 
described in Appendix A[ 

Our anelastic LES code uses the specific entropy s and the three-dimensional 
velocity field u as prognostic variables Pauluis, 2008 . We use a second-order central 
difference spatial discretization scheme with strongly stability preserving Runge-Kutta 
timestepping Spiteri and Ruuth, 2002 . The time step is dynamically adjusted to 


ensure a Courant number near 0.3. Because LES merely resolves the most energetic 
scales of the flow, subgrid-scale (SGS) motions must also be modeled, and we do so 
by applying a constant eddy diffusivity of iz = 1.2 m^ s“^ throughout the domain. We 
choose a constant diffusivity to avoid the nonlinearities that appear in the computation 
of the eddy viscosity by more advanced SGS schemes such as the Smagorinsky-Lilly 
closure Smagorinsky, 1963 Lilly 1962 , which would need to be linearized in a QL 
simulation; constant diffusion as an SGS closure allows a more direct comparison of 
fully nonlinear and QL simulations, notwithstanding that it is an inferior SGS closure. 
The domain extends 6.4 km x 6.4 km in the horizontal and 3.75 km in the vertical, with 
a horizontal and vertical resolution of 25 m. Horizontal boundary conditions are doubly 
periodic. At the upper boundary, flow fields are linearly relaxed over a 800 m deep 
layer toward the horizontal mean flow, which is almost motionless (horizontal mean 
velocities are of the order 0.1 ms“^). The relaxation coefficient varies from r = 0 at 
the bottom of the layer to r = (100 s)“^ at the top. 

The initial state is stably stratified with a potential temperature 9 that increases 
linearly with height, at a rate of 2 Kkm“^. Here, the potential temperature is 
related to the specific entropy by 0 = rexp[(s — s)/cp], where Cp is the specific 
heat at constant pressure, T is a standard temperature, and s is a standard specific 
entropy at the standard temperature T and standard pressure p = 1000 hPa0 The 
initial state is destabilized by imposing a constant sensible heat (enthalpy) flux of 
70.46 W m“^ at the lower boundary. Normally distributed random fluctuations of the 
potential temperature with amplitude 0.1 K in the lowest 200 m break the horizontal 
homogeneity of the initial state and allow the generation of turbulent motions. The 
initial flow is uniform, horizontal, and has a speed of 0.01 ms“^. Together with the 
drag at the lower boundary, this allows for turbulent momentum fluxes to develop 
We ran simulations for up to 12 simulated hours, over which a dry 


Soares 2004 


convective boundary layer forms and grows as a result of the heating at the bottom. 
Because of the statistical horizontal homogeneity of the flow, we use the horizontal 

^ Using the ideal-gas law, it can be verified that 9 is the usual potential temperature with reference 
pressure p. 9 = T{p/po)^^‘^p, with specific gas constant R. However, this potential temperature is 
evaluated at the anelastic reference-state pressure po(z) rather than at the in situ pressure p, as is 


required for thermodynamic consistency of the anelastic approximation [Pauluis 


2008 
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average to define mean fields and eddies, as in Eq. The QL truncation (22), 
here with the state vector '5' = {u,v,w, , is implemented by removing at every 

time step the nonlinear eddy-eddy interactions from the tendencies of all prognostic 
variables!^ [Herring 1963 similarly used the QL approximation to study thermal 
convection between two parallel horizontal plates. 


3.2. Results 


Figure [2 compares vertical and horizontal cross sections of the vertical velocity field in 
the fully nonlinear and in the QL simulations. The vertical cross sections (Fig. [^,b) 
show that upward motion mainly occurs in vertically coherent updrafts, as is well 


known [e.g., Schmidt and Schumann 1989 . In the QL simulation, updrafts are more 
coherent because small-scale structures are missing while the larger scales are well 
represented. The horizontal cross section in the fully nonlinear simulation reveals 
that the updrafts are organized into polygonal convective cells (Fig. &)■ Such 
cellular patterns are well known to arise near the onset of thermal instability of a 


fluid heated from below Chandrasekhar, 1961 chapter 2]. The QL simulation does 
not capture these horizontal correlation patterns (Fig. Si), probably because eddy- 
eddy interactions in the horizontal plane play a role in generating the smaller scales 
of the cellular patterns. The vertical velocity fluctuations w' are distributed more 
symmetrically around zero in the QL simulation: updrafts and downdrafts are of 
similar size and strength (Fig. [^, d). This stands in contrast to the fully nonlinear 
simulations, in which updrafts are faster and narrower and downdrafts slower and 
broader. 

Figurej^ shows the evolution of the vertical profile of the potential temperature in 
the full and QL simulations. Because of the constant heating at the lower boundary 
and the absence of any thermal relaxation in the fluid interior, the boundary layer 
continually deepens and does not reach a statistically steady state (Fig. [^). The 
boundary layer is well mixed with homogenized potential temperature below its top. 
This indicates that the BL is close to the critical state for thermal instability. The 
QL simulation captures quite accurately the growth rate of the boundary layer and 
the mixing of potential temperature below its top. This is because the growth of the 
boundary layer mainly arises through interactions of fluctuations with the mean flow, 
which are retained in the QL simulation. 

Above the top of the boundary layer, defined as the altitude below which 
potential temperature is well mixed, the QL and the fully nonlinear simulations exhibit 
important differences. In the QL simulation, the mean potential temperature profile 
is identical to the initial prohle. By contrast, the fully nonlinear simulation shows a 

® The QL truncation \22\ is valid for the Boussinesq approximation. The anelastic approximation 
requires us to use an average weighted by the background density, as explained in [Appendix A| But 
because averages are here performed on horizontal surfaces with constant background density, the 
QL approximation for the anelastic system is also given by ||22|l. 
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Figure 1. Instantaneous cross sections of the vertical velocity after 4 hours. 
Vertical cross sections in (a) fully nonlinear simulation and (b) QL simulation. 
Horizontal cross sections at 250 m altitude in (c) fully nonlinear simulation and 
(d) QL simulation. Note the different color scales for the fully nonlinear and QL 
simulation. 


layer of strong stability associated with convective overshoots of thermal updrafts into 
the free atmosphere and the downward entrainment of warmer free-atmospheric air 
into the boundary layer [e.g., de Roode et al. 2004 . These overshoots are missing in 
the QL simulation, as they are in first-order diffusive closure schemes for convective 
boundary layers [e.g., Stull 1988 . The difference at the top of the boundary layer is 
emphasized in the profile of temperature variance (Fig. If) , which shows a strong peak 
above the top of the boundary layer for the fully nonlinear simulation but vanishing 
variance throughout the mixed layer and aloft for the QL simulation. Near the surface, 
the mean potential temperature is reduced in the QL simulation, which reflects a 
more efficient upward transport of heat into the interior of the boundary layer. This 
probably can be ascribed to the more intense vertical velocities and the strengthened 
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Figure 2. (a) Vertical profiles of potential temperature at indicated times in the 
fully nonlinear (solid lines) and QL (dashed lines) simulations, (b) Corresponding 
vertical profiles of the potential temperature variance 

a b 




Figure 3. Turbulence kinetic energy budget of the fully nonlinear (a) and QL 
simulations (b). The different terms denote the TKE advection TKE transport 
by eddies T, shear production tS, buoyancy production pressure correlation 
term P, dissipation to subgrid-scales T> and rate of change of TKE 7^. Also 
shown are the residuals, which are of the same order as the SGS diffusion term in 
the fully nonlinear simulation [cf. [Heinze et al.| |2015] but are smaller in the QL 
simulation, likely because the latter is smoother, thus limiting numerical diffusion. 


correlations between buoyancy and vertical velocity fluctuations in the more coherent 
updrafts of the QL flow. 

An inability of diffusive closures and the QL simulation to represent the turbulent 
transport of turbulence kinetic energy (TKE) lies behind the failure of the QL 
simulation and diffusive closures to represent convective overshoots. We take TKE 
to be the kinetic energy of the resolved scales of the LES: e = 0.5 
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Its mean budget is derived from the momentum equations and reads 


dte = —w dzC - dz (pow'e) — w'u'dzU + w'v'dzV + w'w'dzw] 

Po ■' -„-^ 

-4 '---' ,c 


T 


+ w'b' —^dz iw'p') +e. (23) 

V 


The right-hand side contains all terms that generate, destroy, or redistribute TKE: 
advection by the mean flow A and the eddies T, shear production S, buoyancy 
production B {b' = gO 'denotes buoyancy fluctuations), the pressure correlation 
term V and dissipation e. The dissipation e denotes the energy flux from the resolved 
scales to sub-grid scales, which in our case is parameterized as viscous dissipation of 
the resolved fields with constant eddy viscosity v. All but the turbulent transport 
term T = — IjpQdz (pow'e) are of second order in fluctuations and hence are retained 
in a QL truncation. 

The TKE budgets for the fully nonlinear and QL simulations are shown in 
Figure The dominance of the buoyancy production term relative to the shear 
production term in both plots indicates that the flow is thermally driven. In the fully 
nonlinear simulation, buoyancy production is positive throughout the well-mixed part 
of the boundary layer (i.e., upward buoyancy flux), zero at the top of the boundary 
layer, and negative aloft (i.e., downward buoyancy flux). The negative flux is related 
to the overshooting thermals, which trigger a downward entrainment flux of free 
atmospheric air into the boundary layer. This negative flux is missing in the QL 
simulation, in which the buoyancy flux is zero above the top of the boundary layer. 
However, the buoyancy flux is well captured in the interior of the boundary layer. 

The triple correlation transport term T represents the vertical transport of TKE 
by eddies. In the fully nonlinear simulation, eddies transport TKE from lower levels 
with high TKE to higher levels with low TKE. This transport seems to be crucial for 
the growth of TKE in the upper part of the boundary layer and especially across the 
top of the boundary layer. The neglect of this term in the QL dynamics is responsible 
for the missing overshoots of thermals across the boundary layer top and the missing 
associated negative buoyancy flux. The transport term T can be evaluated in the 
QL simulations and actually is nonzero. However, it decouples from the second-order 
dynamics and so does not affect the TKE budget. 


3.3. Implications 

These results show that a QL simulation can capture important aspects of the 
evolution of a dry convective boundary layer, such as its well-mixed nature and its 
growth over time. They suggest that a corresponding CE2 closure would also capture 
the relevant first-order statistics and, therefore, has promise as a nonlocal second-order 
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closure in climate models. Deficiencies such as the missing convective overshoots at 
the top of the boundary layer may be remedied, for example, by adding a linear, 
diffusive transport term in the prognostic equations of the second-order moments [e.g. 
Mellor, 1973 . The resulting parameterization scheme would solve directly for the 1st¬ 


and 2nd-order statistics in every grid box of the large-scale model. 

The applicability of such a scheme depends on whether the corresponding 
CE2 equations can be solved efficiently—much more efficiently than an explicit QL 
simulation can be run. This will require a simplified representation of horizontal 
covariances, for example, by assuming approximate statistical symmetries such as 
horizontal homogeneity and isotropy of fluctuations. But the important nonlocal 
effects of vertical covariances need to be retained. 

The eddy-eddy interactions neglected in the QL simulations may be even more 


critical for moist convection than they are for the dry convective boundary layer Firl 


and Randall 2015 . Hence, more sophisticated approaches may be needed to expand 


the field of application of CE2 closure schemes to moist boundary layer dynamics. 
Exploring to what degree the CE2 approximation and its extensions can capture the 
dynamics of clouds and moist boundary layers promises to be a fruitful area of study. 


4. Large-scale eddy decay on the rotating sphere 


While the preceding example concerned the applicability of CE2 to atmospheric 
dynamics on scales of meters to kilometers, we now turn to a prototype problem 
for atmospheric dynamics on scales of hundreds to thousands of kilometers. Eddies on 
such large scales are essentially the well-known weather systems. They are generated 
by baroclinic instability and are fundamental for maintaining Earth’s climate because 
they are responsible for the bulk of the atmospheric transport of energy, water vapor, 
and angular momentum. Through these transports, they shape the distribution of 
temperature, precipitation, and winds at the surface Peixoto and Oort 1992 . The 


fundamental balances governing such large-scale eddies are different than those in the 
boundary layer. The Coriolis force due to the planetary rotation and the average 
stable stratification become of primary importance, leading to flows that are more 
two-dimensional in character than boundary-layer flows. A two-dimensional (latitude- 
longitude) model suffices to illustrate some of the issues one encounters if one wants 
to develop a closure for the large-scale dynamics of the atmosphere based on cumulant 
expansions. 


4 . 1 . Barotropic model for the upper troposphere 


Large-scale eddies in Earth’s atmosphere are generated near the surface in 
midlatitudes, propagate upward through the troposphere, and propagate meridionally 


in the upper troposphere Simmons and Hoskins 

1978 

Edmon Jr et al. 

1980 

Held 

and Hoskins 

1985 

Thorncroft et al. 

1993 . Their meridional transport and eventual 
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dissipation by wave breaking in latitude bands away from their generation latitudes is 
what gives rise to their meridional angular momentum transport: Large-scale eddies 
in rapidly rotating atmospheres transport angular momentum from their dissipation 
latitudes into their generation latitudes, that is, in the opposite direction to their 
meridional propagation Kuo 1951 Held 1975, 1999 . This angular momentum 


transport ultimately shapes the strength and distribution of surface winds, with 
easterlies in the tropics, westerlies in midlatitudes, and weak easterlies again in 


polar latitudes [see Schneider 2006, for a review]. To understand the strength and 
distribution of surface winds, it is therefore necessary to understand the meridional 
propagation and dissipation of large-scale eddies, which are concentrated in the upper 
troposphere Ait-Chaalal and Schneider 2015 . The simplest model that captures 


these processes is the barotropic model—a model of a two-dimensional fluid layer on a 


sphere, thought to represent a layer in the upper troposphere [e.g.. Held and Phillips 


1987 


Equations of motion The equations governing barotropic flow can be derived from 
the Boussinesq equations 0. Consider a Boussinesq flow on a sphere of radius a 
rotating at constant spin angular velocity fJ. We further assume that the flow is two- 
dimensional on the sphere: u • = 0, with being the radial coordinate. In that 


case, the momentum and continuity equations (la lb) reduce to 


dv 

— —I- V • Vv -I- 211 X V = — Vd) + Ev 
at 

V-v = 0. 


(24a) 

(24b) 


The horizontal components of the wind and of the forcing are denoted v and Jv- 
Because the flow v is two-dimensional on the sphere, only the local vertical component 


SI sin (latitude 4>) of H yields non-zero terms in (24a). Taking the curl of the 
momentum equation and projecting it onto the radial direction e^. yields the two- 
dimensional barotropic vorticity equation, 


dq 


V • Vg = (V X Ev) ■ Gr- (25) 

The flow is entirely described by this equation for the absolute vorticity 

<7 = / + C, (26) 

where the Coriolis parameter /((/>) = 2S1 sincj) represents the vorticity of solid body 
rotation, and C = (V x v) • is the relative vorticity in the radial direction e^, 
relative to the rotating reference frame. The advection term v • Vq contains the 
advection of planetary vorticity v • V/ = fSv, with /3 = a~^d^f = 2Sla“^ coscj), which 
arises from the curl of the Coriolis force. This term is commonly referred to as the 


/3-term. The vorticity equation (25) contains the entire dynamics of the flow because 
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an incompressible two-dimensional flow is described by a streamfunction ip, defined 
such that 


V = V X (ipBr), 

c = VV- 


(27a) 

(27b) 


That is, if the relative vorticity (p is known, the advecting velocity (27a) can be 


determined from the streamfunction, which is the solution of a Poisson equation (27b). 


Thus, the equations of motion (25) and (27), supplemented by appropriate boundary 


conditions, specify the dynamics completely. 


The equations of motion (25 )-(27) can be nondimensionalised using the planetary 
radius a as the typical length scale and the length of the day 27rr2“^ as the typical 
time scale. With that nondimensionalisation, the operators V, Vx and become 
operators on the unit sphere, and the angular velocity VI becomes 27r. Throughout the 
rest of the paper, we will use there nondimensionalized quantities, unless otherwise 
specihed. 


Eddy-mean flow decomposition We consider situations in which the boundary 
conditions of the problem are statistically symmetric under rotations around the 
planet’s spin axis, so that the flow statistics (but not the instantaneous flow itself) can 
be expected to be axisymmetric. A zonal average (•) then suggests itself. Decomposing 
flow fields in the nondimensional barotropic vorticity equation into zonal means (•) 
and eddies (•)' = (•) — (•) yields the mean and eddy vorticity equations. 


(K 

dt 


^ = -UiP',C'), 


d 


d 


dt cos (p dX 


c = - 






ac 

d(p 


(28a) 

(28b) 


Here (3 = 2flcos(p, with H = 27r, is the nondimensional meridional derivative of the 
Coriolis parameter, and the Jacobian 




1 / dip dtp dip dp \ 

cos (p \dX d(p d(p dX J 


(29) 


represents the advection on the unit sphere of vorticity () by the zonal (u) and 
meridional (n) flow implied by the streamfunction ip: 


dip 


1 dip 
cos (p dX' 


(30) 


The streamfunction-vorticity relations are 


C = 
C' = 


(31a) 

(31b) 
















20 


with V„ likewise defined on the unit sphere. 

Equation (28) is essentially (19) for the 2D barotropic vorticity equation. 
The mean flow evolves in time due to vorticity fluxes — The eddy 

vorticity budget involves shear by the mean flow uj cos4)d\(^', eddy-eddy interactions 
Jn{ip'X') ~ Jn{fp'X') ; the /?-term (3v' and the advection of mean-flow vorticity 
dXv'. 


Cumulants The first cumulant is the mean vorticity t), and the second cumulant 
is the vorticity equal-time two-point correlation: 


CXX) = 


c(ri,r2,t) = = C(0i>2, Ai - 


(32a) 

(32b) 


The first cumulant depends on the latitude (j), and the second cumulant depends on 
the latitudes 4 >i and 4>2 and, because of the statistical axisymmetry of the equations, 
on the longitude difference Ai — A 2 [Marston et al. 2008 . Because the flow is entirely 
determined by the scalar q (or Q, all other correlations are determined by the scalar 


cumulant c. Hence, moments like ^'(ri,f) 
calculated from ( and c [e.g., 


u'(r 2 ,t) or u'(ri, f) (g) u'(r 2 , t) can be 


Marston et al. 


2008 


Srinivasan and Young 


2012 


Marston et al., 2014 


The CE2 equations for this problem are given in [Marston et al. 2008, 2014 


They are of the form (I 81 ) and (18:), with vorticity fluxes appearing as the essential 


eddy terms. The equivalent QL system is (28) where the eddy-eddy interactions are 
neglected. 


Numerical implementation We simulate a barotropic flow on a sphere, specified by 
the equations of motion ([^ and ([^ on a spherical geodesic grid [Heikes and Rand^ 


1995a|b Qi and Marston , 2014 with 163,842 cells. To remove enstrophy that cascades 
to the smallest scales, hyperviscous dissipation j^(V^ -I- 2)V®() is included, where the 
term (V^ -I- 2) ensures that angular momentum is conserved. The hyperviscosity 
coefficient ly is chosen such that the smallest resolved mode decays with an e-folding 
time of 5. The vorticity is stepped forward in time by a second-order leapfrog scheme 
using the Robert-Asselin-Williams filter Williams, 2009 . The time step is fixed at 
At = 0.01. 


Explicit time integration of the cumulant equations is carried out in spectral 
space using a 4th-order Runge-Kutta algorithm with an adaptive time step. We 
adopt the spectral truncation 0 < £ < L on spherical wavenumber i, with the zonal 
wavenumbers restricted to \m\ < min{£, M}. We choose spectral cutoffs L = 60 and 
M = 30. Hyperviscosity is adjusted to ensure the same e-folding time at the maximum 
wavenumber £ = L as on the smallest resolved spatial scales on the geodesic grid. 

To verify that the spectral cumulant simulation has sufficient resolution and can 
be compared to the geodesic grid model, a simulation of the fluid is also performed 
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in a pure spectral calculation with the same numerical methods and resolution as for 
the cumulant equations. The agreement between the spectral and geodesic models is 
excellent. QL simulations are performed in spectral space by removing the triads that 
correspond to the interaction of two eddies, each with non-zero zonal wavenumber. 

A program that implements fully nonlinear simulations on the spherical geodesic 
grid, and the nonlinear, QL, and CE2 simulations in spectral space, is freely available]^ 
More details about the simulations and the cumulant expansions can be found in 


Marston et al. 


2014 


4-2. Eddy lifecycle simulations 


Setup To illustrate situations when CE2 and QL approaches succeed or fail at 
capturing barotropic flow dynamics, we simulate the evolution of an initial zonal flow 
U{(f) with a superimposed initial disturbance (eddy) with vorticity Q{4>, A). The zonal 
flow U and disturbance Q mimic the upper-tropospheric jet stream and disturbances 
that may originate, for example, from lower-tropospheric baroclinic instability. The 
setup is inspired by Held and Phillips 1987 and uses 


U{4)) = A cos (j) — B cos^ <l) + C sin^ </> cos^^ cj) — D cos® cj), (33a) 

Cii4>,X) = Cocos(l)exp[-{(j)-(j)m)^/6'^)]cos{kiX). (33b) 


To mimic Earth’s upper troposphere, we choose 

A = 3.4 X 10-\B = 4.1 X 10"\C' = 4.0, and D = 2.3 x 10"^ (34) 


The corresponding dimensionalized flow on a sphere of Earth’s radius and rotation rate 
resembles the midlatitude jet in the upper troposphere. It has a maximum westward 
velocity of ~33ms“^ at a latitude of ~40° and a maximum eastward velocity of 
^22ms“^ at the equator; its zero is near ~17°. This zonal flow is barotropically 
stable. The eddy vorticity Q decays meridionally away from its maximum absolute 
value Co cos (jjm at latitude (j)m = 7r/4 = 45°, with characteristic meridional decay scale 
(5 = 7r/9. Its zonal wavenumber = 6 is close to the dominant zonal wavenumber 


of transient eddies on Earth [e.g., [Boer and Shepherd 1983 Straus and Ditlevsen 


1999 , which approximately coincides with the baroclinically most unstable zonal 


wavenumber [e.g.. 

Simmons and Hoskins 

1976 

1978 

Thorncroft et al. 

1993 Merlis 

and Schneider 

2009 . 


We let the flow evolve freely without forcing or dissipation, apart from 
hyperviscosity, and analyze the time evolution of the mean flow and the eddies. We 
compare CE2 to the statistics of fully nonlinear simulations for different choices of 
parameters. To identify relevant nondimensional parameters controlling the evolution 
of the flow and the adequacy of CE2 closures, we rescale the mean and eddy 

^ The application “GCM” is available for OS X 10.9 and higher on the Apple Mac App Store at 
URL http://appstore.com/mac/gcm 
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vorticity equations (28), using the relative vorticity 2Ro of the initial zonal flow U. 
Dimensionally, we have Ro = A/(2D), where A « 2 max([/)/a is the typical initial 
zonal-mean flow vorticity. Hence, Ro is a Rossby number measuring the mean flow 
vorticity to the equatorial planetary vorticity 212. We measure the initial maximum 
vorticity Co of the eddies relative to the mean-flow vorticity 2Ro through the amplitude 
parameter e = Co cosC)to/ 2 Ro. The quantities C, C^ ■*/'^ ^od t are then rescaled with 
47 rRo, 47 reRo, 47 rRo, 47 reRo, and ( 47 rRo)“^, respectively. The equations of motion for 
the mean-flow and the eddies become 


cK 

dt 


^=-eV„(V^CC'), 


Ro 


d \ dij} d 
dt cos (/> dcj) dX 


C = -eR0 Jr,{i;',C)-JnWX') 


(35a) 

A , ^ 1 9C 

aA ^ ° cos (I>d4> 

(35b) 


Equation (35) gives the relative amplitude of the different terms if we assume that 


the typical length scales of the mean flow and eddies are of the order of the radius 
of the planet. An immediate simplification that results for small Rossby number Ro 
(the case we will consider) is that the advection of mean-flow vorticity C by meridional 
velocity fluctuations is negligible compared with the /3-term, which is a factor Ro“^ 
larger. 

The nondimensional parameters Ro and e control the vorticity of the mean flow 
and of the eddies and are important for the evolution of the barotropic flow. Eddy- 
mean flow interactions are of order Ro, and eddy-eddy interactions of order eRo, 
provided is of order one. This is the case initially; however, it does not remain 
true over the evolution of the flow, as small scales are generated. The two parameters 
Ro and e can be varied independently in our setup. In what follows, we explore how 
these parameters affect the flow evolution and the adequacy of CE2 and QL approaches 
in capturing it. 


4-3. Results 

Varying eddy amplitudes For a fixed initial mean-flow Rossby number Ro « 0.06 


(corresponding to the Earth-like parameters in equation (34)), we run eddy lifecycle 
experiments for larger-amplitude initial eddies with e Ri 6, and for smaller-amplitude 
initial eddies with e « 2. The expectation is that CE2 and QL approaches are 
more successful for the smaller-amplitude eddies, for which the nonlinear eddy-eddy 
interactions (of order eRo) are weaker, and this is indeed borne out in the simulations. 
It is instructive to see in what way they fail to capture aspects of the flow evolution 
for the larger-amplitude eddies. 

For the larger-amplitude eddies, Fig.|^shows the relative vorticity C at 5 instants 
during the evolution of the flow. It is evident that the initial disturbance quickly 
becomes nonlinear and develops drawn-out filaments on the equatorward flank of the 
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zonal jet. The filaments roll-up anticyclonically within cats’ eyes structures (marked 
by X’s in Fig. that continue to have the initial zonal wavenumber ki = 6. Such 
cats’ eyes are characteristic of Rossby waves that break in “surf zones” around their 


critical layerq^l [Stewartson 

1977 

Warn and Warn 

1978 

McIntyre and Palmer 

1983 

Killworth and McIntyre 

1985 

Held and Phillips 

1987 . 



T = 4.0 


Figure 4. Evolution of relative vorticity in the fully nonlinear simulation for the 
larger-amplitude eddies (e = 6) for an Earth-like setting (Ro = 0.06). The relative 
vorticity maps show the formation of cats’ eyes, with vorticity filaments rolling 
up within them. White X’s mark the centres of some cats’ eyes. Time scales are 
non-dimensionalized with the length of the day 27rQ,~^. 

As the eddies break and eventually dissipate in the surf zone, they are absorbed 
by the mean flow, and their kinetic energy decays. The total eddy kinetic Ek = 
0.5 tK cos where ck = 0.5(u'^ becomes very small at large times (> 30, 

see Fig. [5f)- At those times, most of the initial kinetic energy has been transferred 
to the mean zonal flow. The local zonal kinetic energy ez = 0.5u^ (u = of the 

mean zonal flow increases in the core of the midlatitude jet, roughly in proportion 
to the decrease of the eddy kinetic energy ck (Fig. &)■ Total energy Ek + Ez, 

^ Fig. shows some spurious oscillations in the critical layer region. This is the consequence of a 
too weak hyperdiffusion. With hyperdiffusion strong enough to remove the ripples, we observed have 
noticable hyperdiffusive wave absorption over the time scales considered. This obscures the wave 
absorption due to eddy—eddy interactions and blurs the differences between fully nonlinear and CE2 
simulations. Removing the spurious ripples while showing sharp differences between fully nonlinear 
and CE2 simulations would have required a much higher resolution, or perhaps a different numerical 
advection scheme. 
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Figure 5. Evolution of kinetic energy and eddy momentum flux convergence 
(EMFC) in the fully nonlinear simulation (left column) and in a direct CE2 
calculation of the statistics (right column) for an Earth-like setting (Ro = 0.06) 
and for the larger-amplitude eddies (e = 6). (a, d) Eddy kinetic energy e^. (b, 
e) Zonal kinetic energy ez- (c, f) EMFC. Vertical blue dashed lines indicate at 
which times relative vorticity snapshots are shown in Fig. Time scales and 
length scales are non-dimensionalized with the length of the day 27rn~^ and the 
planet radius a. 


with Ez = 0.5 J^ez cos (j)d(j), is approximately conserved in the model, up to very 
small losses (^0.2% of the total after a time of 50) through subgrid-scale dissipation. 
That is, although dissipation at small scales in the surf zone is essential to generate 
irreversibility, the kinetic energy loss associated with it is small compared to the 
transfer to the mean flow. 

The transfer of Ek to Ez implies acceleration of the mean zonal jet. This 
acceleration occurs through transfer of momentum from the eddies to the mean flow, 
as can be seen from the nondimensionalized zonally averaged momentum equation 
(24a) in the inviscid limit (Jv = 0): 


at cos 0 (70 


(36) 


Acceleration of the mean zonal flow occurs where eddy momentum fluxes u'v' cos (j) 
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converge. Multiplying the mean zonal momentum equation (36) by u and integrating 
by parts yields the equation for the zonal kinetic energy Ez, 


dt 


,Ez = 


COS^ (jju'v' 


d 


dcj) ycoBcj) 


d(j) = 


(37) 


where the right-hand side is obtained from a corresponding integral of the eddy 
momentum equations. This shows that transfer of kinetic energy from the eddies 
to the mean flow occurs through eddy momentum fluxes that are up the gradient of 
angular velocity ucos“^ 4>. The acceleration of the mean zonal jet at its core (Fig. &) 
thus is associated with eddy momentum flux convergence d^{u'v' coscj)) (EMFC, see 
Fig.js}?). 

However, the eddy kinetic energy does not decay monotonically. Instead, it 
exhibits damped oscillations during which eddy momentum fluxes cause zonal angular 
momentum to slosh back and forth meridionally within the jet (Fig. §). The 
alternating poleward and equatorward momentum fluxes (with decreasing amplitude) 
on the equatorward flank of the jet are result of nonlinear processes within the 
surf zone. These processes have been described in an idealized analytical model of 
Rossby-wave breaking in critical layers, the Stewartson-Warn-Warn (SWW) theory 
(Stewartson 1977[ Warn and Warn 1978 see also Kill worth and Mclntyre[ 19851. The 
oscillation on the poleward flank of the jet may be more linear, originating from the 
reflection of Rossby wavepackets from reflecting levels that arise because /3 decreases 
with latitude [e.g., 


Lorenz 2014 


Figure]^ (right column) shows the kinetic energies and EMFC obtained from a 
direct calculation of these statistics with CE2. CE2 captures the oscillation of kinetic 
energy between eddies and the mean zonal flow, with a period similar to the fully 
nonlinear simulation (Fig. [5}i, b). However, the oscillations are too weakly damped; 
large eddy kinetic energies ex persist for a long time. The eddy absorption in the surf 
zone is not adequately captured by CE2 because CE2 cannot resolve the generation 
of small scales in the surf zone through eddy-eddy interactions. Consistently, 
unrealistically strong oscillations persist in the EMFC under CE2 (Fig.[^). How these 
oscillations arise from the perspective of wave mechanics, and why CE2 cannot capture 
the wave absorption in this case, is illustrated in [Appendix B| in a QL simulation 
that corresponds to the CE2 calculations shown here. The phenomenology of such 


oscillations has been described analytically by Haynes and McIntyre 1987 in the 
context of a QL truncation of the SWW model. 

Eddy-eddy interactions are weaker and CE2 is more successful in capturing the 
flow dynamics when the amplitude of the initial perturbation is decreased by a factor 
3 (e « 2). Cats’ eyes with rolling-up vorticity filaments do not develop in the fully 
nonlinear simulation (Fig. [^. Instead, eddies are sheared by the mean flow, which 
transfers eddy kinetic energy ex to the mean flow through the Orr mechanism, which 
is only weakly nonlinear because it involves the interaction of disturbances with the 
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Figure 6. Evolution of relative vorticity in the fully nonlinear simulation for the 
smaller-amplitude eddies (e = 2) for an Earth-like setting (Ro = 0.06). Time 
scales are non-dimensionalized with the length of the day 2'kQ,~^ . 


slowly varying mean flow [e.g. Thomson, 1887 Orr 1907 Farrell, 1987 Bouchet et al. 
2013 . The transfer of eddy kinetic energy to the mean flow occurs over time scales of 


a couple days, corresponding to the shear time scale of the mean zonal flow (Fig. [^. 
The damped oscillatory behavior seen in the larger-amplitude simulation disappears. 
Because eddy absorption results from the mean flow shearing the eddies—an eddy- 
mean flow interaction that is captured by CE2—statistics calculated directly with CE2 
come in very close agreement with those from the fully nonlinear simulation (Fig. 
right column). As in the nonlinear case, eddy absorption occurs through the formation 
of small-scale vorticity filaments. But instead of rolling up inside cat’s eyes, here they 
stretch around the planet. 

It is worth noting that eddies are also sheared equatorward of the cats’ eyes in 
the larger-amplitude simulation (Fig. |^. Hence, weakly nonlinear eddy absorption 
also occurs in this simulation, but it is not the dominant effect responsible for eddy 
absorption. 


Relative amplitude of terms in the potential vorticity budget For the larger-amplitude 
experiment, the evident importance of the development of small scales through eddy- 
eddy interactions seems consistent with the order of magnitude of the terms in the 
vorticity equations (35). The eddy-eddy interactions appear of order eRo « 0.4, 


compared with the ,d-term which is responsible for Rossby wave dynamics and is of 
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Figure 7. Evolution of kinetic energy and eddy momentum flux convergence 
(EMFC) in the fully nonlinear simulation (left column) and in a direct CE2 
calculation of the statistics (right column) for an Earth-like setting (Ro = 0.06) 
and for the larger-amplitude eddies (e = 2). (a, d) Eddy kinetic energy e/f. (c, 
f) EMFC. Vertical blue dashed lines indicate the times corresponding to relative 
vorticity snapshots on Fig.|^ Time scales are non-dimensionalized with the length 
of the day 2nfl~^ and length scales with the planet radius a. 


order 1. Hence, the eddy-eddy interactions are not negligible compared with Rossby 
wave dynamics. By contrast, the interactions of the disturbance with the mean flow 
shear are of order Ro « 0.06 and hence are much weaker. 

For the smaller-amplitude experiment, the dimensional analysis of the vorticity 
equations (35) suggests that the eddy-eddy interactions now are of order eRo < 0.2 
compared with the /3-term. Hence, the eddy-eddy interactions become close to 
being negligible compared with Rossby wave dynamics, consistent with the simulation 
results. However, the dimensional analysis also suggests that interactions of the 
disturbance with the mean flow shear still are of order Ro « 0.06 and hence are 
weaker still, albeit of the same order of magnitude as the eddy-eddy interactions. Yet 
the eddy-eddy interactions are inhibited in the fully nonlinear simulation, whereas the 
shear interactions dominate the eddy absorption, illustrating the limits of dimensional 
analysis in this nonlinear problem. 

Equation (35) is useful to determine which parameter to vary. Nevertheless, one 
has to be careful when interpreting the relative amplitude of the different terms. A 
small term can be fundamental for the dynamics. For example, the SWW theory 
has shown that eddy-eddy interactions, albeit weak, can determine at leading order 
momentum fluxes because they dominate the vorticity budget in a thin critical layer. 
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where linear theory breaks down. This remains true in the limit where the relative 
amplitude of the eddy-eddy interactions goes to zero. Moreover, the relative amplitude 
of the terms in the vorticity budget evolves with time as small scales develop. 


Varying Rossby numbers To illustrate how variations of the mean-flow Rossby 
number affect the evolution of disturbances, we use larger-amplitude eddies (e « 6) and 
decrease the Rossby number Ro from 0.06 to 0.02. Dimensionally, this is equivalent 
to weakening the initial flow while keeping the rotation rate of the planet constant, 
or to increasing the rotation rate while keeping the initial flow constant. Based on 


the dimensional analysis of the vorticity equations (35), this reduction of the Rossby 


number should decrease the relative magnitude both of eddy-eddy interactions and of 
shearing of eddies by the mean flow relative to the /3-term, maintaining the relative 
amplitude of the eddy-eddy interactions to the /3-term. 

Figure [^compares the time evolution of the eddy kinetic energy ck in the fully 
nonlinear and the CE2 simulations. As we have already seen, eddy absorption at 
Ro = 0.06 is not captured by CE2 because it is strongly nonlinear. However, at the 
smaller Rossby number Ro = 0.02, it is faithfully captured by CE2. Indeed, eddy 
absorption appears more weakly nonlinear for Ro = 0.02, and dominated by mean- 
flow shearing, as is evident on the relative vorticity maps (Fig. [^. Vorticity maps 
for Ro = 0.04 and Ro = 0.03 show the transition from weakly to strongly nonlinear 
absorption: the meridional extent of the nonlinear surf zone decreases with increasing 
rotation, while shearing effects become more important (Fig. [^. Because CE2 can 
capture the weakly nonlinear shearing interactions but not the strongly nonlinear 
eddy-eddy interactions in the surf zone, it becomes gradually more adequate as the 
rotation rate increases (Fig. |^. This occurs although the orders of magnitude of the 
relevant terms suggested by the dimensional analysis decrease by equal factors of order 
Ro as the mean-flow Rossby number decreases. It appears that what is important here 
is that the magnitude of the advection of absolute vorticity, and hence of linear Rossby 
wave dynamics, increases relative to both of these terms. For Ro = 0.02, shear explains 


eddy decay and jet acceleration, even though the nondimensionalization (35) suggests 
shear should be much smaller than eddy-eddy interactions. 


Higher-order closures CE2 fails to capture eddy absorption when eddy-eddy 
interactions are important for the dynamics. We tested whether a higher-order closure 


(CE3*) captures eddy absorption more faithfully. CE3* is described in Marston et al. 


[2014| . It truncates the cumulant equations at third order, ensuring realizability by 
projecting out modes with (unphysical) negative energies. 

Because CE3* is computationally 


The results are summarized in Fig. 10 


expensive, the resolution has been reduced to M = 40 and L = 20. For simplicity we 


turn off eddy damping (t = oo, see Marston et al. 2014 for definitions). The full and 
CE2 simulations in Fig. are run at this lower resolution and are consistent with 
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Figure 8. Evolution of eddy kinetic energies ex in the fully nonlinear simulation 
(left column) and in a direct CE2 calculation of (right column) for the larger- 
amplitude eddies (e = 6), with Rossby number decreasing from Ro = 0.06 to 0.02. 
Time scales are non-dimensionalized with the length of the day 2'kQ,~^ and length 
scales with the planet radius a. 


the higher-resolution runs (Fig. [^; however, they do exhibit a faster eddy damping 
because of the stronger diffusion. It appears that CE3* captures the eddy absorption 
very accurately. We also tested other closures that take into account third-order terms 


Marston et al. 2014 ; they bring similar improvements. 
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T=1.2 T = 3.0 T = 7.5 T=17.5 



Figure 9. Evolution of relative vorticity in the fully nonlinear simulation for the 
larger-amplitude eddies (e = 6), with Rossby number decreasing from Ro = 0.06 
to 0.02. For comparison, the first row reproduces the vorticity maps already shown 
in Fig.|^ Time scales are non-dimensionalized with the length of the day 27rQ“^. 

4.4- Implications 

The results show that direct CE2 calculation of barotropic flow statistics representative 
of the upper troposphere can succeed in circumstances when the dominant nonlinear 
interaction is between eddies and the mean flow, for example, by shearing. They 
fail when strongly nonlinear eddy-eddy interactions become important in surf zones 
around critical layers, where the roll-up of vorticity filaments leads to the generation 
of small scales. This is a process that cannot be captured in CE2. However, higher- 
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Full 
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Figure 10. Evolution of eddy kinetic energies ex (top) and eddy momentum 
fluxes EMFC (bottom) in the fully nonlinear simulation (left column), in a direct 
CE2 calculation (middle), and in a CE3* calculation (right), all for the larger- 
amplitude eddies (e = 6) and an Earth-like setting (Ro = 0.06). Compared with 
the simulations in previous figures, the resolution is reduced to M = 40 and 
L = 20. 


order closures, which take some effects of third cumulants on second cumulants into 
account, can perform better in such cases—at the price of increased conceptual and 
computational complexity. 

Weakly nonlinear eddy-mean flow interactions seem to be favored over strongly 
nonlinear eddy-eddy interactions when the eddy vorticity is small enough compared 
with the planetary vorticity. The transition from weakly to strongly nonlinear 
interactions predominating occurs above a critical value of eddy amplitude e that 
is a decreasing function of the Rossby number Ro. The parameter e need not be small 
for absorption through weakly nonlinear shearing to occur. For example, for small Ro, 
weakly nonlinear eddy absorption through shearing seems to be favored even when a 
large value of e suggests that nonlinear eddy-eddy interactions should be larger than 
the mean-flow shearing of eddies. 

When strongly nonlinear eddy-eddy interactions are favored, high eddy kinetic 
energies develop in the QL/CE2 approximation because momentum sloshes back and 
forth meridionally within the jet, without sufficient absorption. Lifecycle experiments 
carried out with a QL GCM have shown that this phenomenology is relevant to the 
upper troposphere in an Earth-like setting, highlighting the relevance of the simplified 
barotropic model: Earth’s upper troposphere appears to be in a regime in which 
nonlinear eddy-eddy interactions in surf zones are important for the structure of the 


momentum fluxes Ait-Chaalal and Schneider 2015 . 
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5. Conclusions 


Atmospheric flows are highly anisotropic and inhomogeneous, with rich spatial 
structures. Turbulent closures that respect the anisotropy and inhomogeneity 


may enable the direct statistical simulation of Earth’s atmosphere Marston 2012 


Expansion of statistics in equal-time cumulants yields equations of motion for the 
statistics that can already provide useful closures at second order (CE2), because the 
mean flows and interactions of perturbations with them are strong [Herring 1963 


[O’ Gorman and Schneider] [2007[ [Marston et al.[ |2008[ [Tobias and Marston[ |2013[ 
Marston et al. 2014 . CE2 solves the first two cumulant (central moment) equations 


of the QL approximation, in which interactions between mean flows and fluctuations 
around it are retained, while nonlinear eddy-eddy interactions are neglected. In 
section[^ we formulated CE2 for the Boussinesq equations by introducing a condensed 
tensorial notation. The case of the anelastic equations is presented in [Appendix A 
and involves the use of density weighted averages. We tested the relevance of CE2 
to two distinct atmospheric flows involving different length and time scales and force 
balances: turbulent convection in the atmospheric boundary layer (section]^, and 
weak two-dimensional turbulence representative of the upper troposphere (section]^. 

Convection in the atmospheric boundary layer links large-scale atmospheric 
dynamics aloft to the surface underneath, mediating the exchanges of momentum, 
energy, and water between the surface and the free troposphere. Motions in boundary 
layers and in clouds that have their roots in them have dynamical scales of meters, 
meaning that they need to be parameterized in GCMs. Current parameterization 
schemes have numerous shortcomings; our inability to represent cloud and boundary 
layer dynamics adequately in climate models is the largest source of uncertainty 
in climate change projections. Because cumulant expansions capture interactions 
between fluctuations (e.g., thermals) and mean fields and take non-local correlations 
of fluctuations into account, without requiring the introduction of tunable parameters 
that proliferate in current parameterization schemes, they may offer a way to achieve 
more physically consistent parameterizations. We presented encouraging initial 
results, showing that a QL large-eddy simulation of a dry convective boundary 
layer captures the first-order statistics (e.g., mean boundary layer growth) of a 
corresponding fully nonlinear LES. However, it does not capture second-order statistics 
adequately. More work is required to investigate to what degree these results 
hold generally, in broader classes of boundary layer flows and in the presence of 
moisture effects and clouds, and how the QL results can be improved by including 
representations of higher-order effects, such as the turbulent transport of kinetic 
energy. 

The potential for development of parameterisation schemes based on GE2 is a 
promising direction for future research. But it requires overcoming both technical and 
theoretical challenges. On the technical level, fast numerical methods are required for 






























33 


the method to be competitive with other subgrid schemes. This could be achieved by 
dimensional reduction to capture only the most important non-local correlations in 
the second cumulant. At the theoretical level, it is not clear to what extent CE2 and 
possible extensions will be able to describe moist convection, or what effects vertical 
shear will have on its accuracy. It appears necessary to include some effects of eddy- 
eddy interactions, such as those captured by the third order CE3* approximation 


Marston et al. 


2014 


At the planetary scale, how and where eddies in the upper troposphere dissipate 
controls the strength and direction of momentum fluxes and thus climatic features 
such as surface winds. A one-layer barotropic model that mimics the behavior of the 
upper troposphere illustrates different mechanisms through which eddy absorption 
by the mean flow can occur. CE2 describes eddy absorption well when it occurs 
through shearing of eddies by the mean flow. This happens when the vorticity that 
characterizes the eddies is small compared with the planetary vorticity (planetary 
rotation rate). When the eddy vorticity is large, CE2 is not adequate because eddy 
absorption results from the formation of small scales that form through eddy-eddy 
interactions in critical layers. A comprehensive theory that describes these weakly 
and strongly nonlinear absorptive regimes is lacking. 

Our results suggest that, in general, higher-order closures are required for an 
accurate direct statistical simulation of large-scale and smaller-scale atmospheric flows. 
We have tested a few of them in the large-scale context and found improvement 
compared with CE2. Nevertheless, going beyond CE2 raises a number of questions. 
Higher-order closures are several orders of magnitude slower than CE2; currently, they 
are much slower than direct simulation of the flows, while CE2 can be faster that direct 
simulations. Dimensional reduction may be able to help here. More generally, once 
eddy-eddy interactions are taken into account, the whole hierarchy of cumulants is 
active and not completely described by any finite closure. Realizability of closures 
becomes an issue [e.g., Orszag, 1970, 1973 Marston et al. 2014 , and it is known that 


intermittency cannot be adequately described in this way Frisch 1995 Lesieur, 2008 


But the existence of anisotropic and inhomogeneous mean flows provides a starting 
point for a systematic exploration of statistical closures. 
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Appendix A. Cumulant expansion of anelastic flow 

Because the flow is non-divergent in the Boussinesq equations, it is possible to directly 
use any averaging operator satisfying the properties 0. However, for more general 
flows, in which the density may vary, the requirements on the averaging operator need 
to be modified so that second-order equations for fluctuations consistent with the 
conservation laws are obtained. Because density appears as a weight in all integrals of 
conserved quantities over the flow domain, the density weighted average (•) = {p-)/p, 
with (•) satisfying the Reynolds properties Q, leads to a decomposition of the flow 
that is amenable to closures such as CE2. The density weighted average satisfies &-c) 
but, importantly, not commutativity with partial differentiation 0 )- Thus, we can 
define an eddy mean flow decomposition 

f = f* + f, (A. la) 

f9* = f*9*+fg, (A. lb) 

-* 

with hats (•) = (•) — (•) denoting fluctuations around the density weighted mean. 
This is sometimes called the Favre decomposition. 

How a varying density affects cumulants and CE closures can be illustrated with 
the anelastic equations, an extension of the Boussinesq equations that allows the 
reference density po = po{z) to vary with height z. In the anelastic approximation, 
density perturbations Sp, the pressure potential $ = 6p/po, and the buoyancy 
b = —g6p/po are defined relative to this height-dependent and hydrostatically balanced 
reference density. The state vector is defined as in 0: but the augmented state 
vector ’4' is now taken to be 


’i' = {u ,v ,w ,b, V$)^ 


(A.2) 


We have to consider covariances involving and not $ because multiplying by 
density does not commute with derivatives. For simplicity, we assume that the linear 
operator £ is homogeneous in the reference density and satisfies poC{'^) = C{po^)- 
It implies that it cannot involve any spatial derivative. The anelastic equations in a 


reference frame rotating with angular velocity fl are Vallis[ 2006 
dpo'^ 


dt 


-I- V • 0 u] = L[po'®'] - PoV$ - 

V • pqu = 0. 


(A.3a) 
(A.3b) 


The Boussinesq equations Q are obtained from the anelastic equations by setting 


po to a constant. The continuity equation (A.3 3) again reduces to a non-divergence 


constraint. Because the reference density po appears in it, we need to use the reference- 

-# - 

density weighted average (•) = (po')/Po to derive the cumulant expansion. 
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The equation for the first cumulant is obtained by averaging (A.3): 


dt 


+ V • [po’®' 0 u*] = -V • [po®' 0 u ] + L[/9 o ^ ] - poV0 


(A.4a) 


V • po u* = 0. (A.4b) 

To define a second central moment or second cumulant, we need a quantity that 


respects the symmetry (111 of the covariance matrix under exchange of the spatial 


coordinates ri and r 2 , that gives the second-order term in the first cumulant equation 


(A.4), and that corresponds to a density weighted average when ri = r 2 . A possible 
choice is 


C+(ri,r2,t) = ^[^ 0 ( 1 - 1 ) +/3o(r2)] ® '®'(r2,t)- 


We also define the two auxiliary covariances corresponding to (12) 


C+(ri,i'2,t) = ^[po(ri) +Po(r2)] $(ri, t)'^(r2, t), 

C“(ri,r2,t) = ^[po(ri)+Po(i’2)]’®'(r2,t)(8)u(ri,t), 


(A.5) 

(A.6a) 
(A.6b) 


and introduce the corresponding quantities, C_, Cj, and C“ with the density 
difference instead of the sum, as in 


C_(ri,r2,t) = ^[po(ri) - ^ 0 ( 1 - 2 )] ’®'(ri,t) 0 4'(r2,t). 
Some algebra then gives the equations of motion for the second cumulant: 
d 


(A.7) 


dt 


C±(ri,r2,t)-I-Vri • [C±(ri,r2,t)(8)u*(ri,t)]-t-Vr2 • [C±(ri, r2, t) 0 u*(r2, t)] = 
- (r 1, r2, t) [vr (r2, t)] ^ (r2, t)C^ (r 1, r2, t) 

- ^C±(ri,r2,t) [V • u*(ri,t) -f V • u*(r2,t)] 

- iCzp(ri,r2,t) [V • u*(ri,t) - V • u*(r2,t)] 

-I- LC±(ri,r2,t) -h C±(ri,r2,t) -f VriC±(ri,r2,t) -t- Vr2C|(r2, ri, t) 

-I- third-order term (A.8a) 


Vr, • [C!;(ri,r2,t) + C“(ri,r2,t)] = 0, 

Vrr • [C!;(ri,r 2 ,t)-C“(ri,r 2 ,t)] = 0. 


(A.8b) 


Truncation of the cumulant expansion at second order (CE2) consists of equations 


(A.4) and (A.8), neglecting the third-order term. Clark and Spitz 1995 have 


derived the evolution equation of the single-time two-point covariance tensor for the 


compressible Navier-Stokes equation, yielding a more general version of (A.8). 
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Appendix B. Eddy absorption in the QL approximation 


To illustrate more clearly why CE2 fails to capture the absorption of larger- 


amplitude eddies, we perform QL simulations of the barotropic vorticity equation (24), 


eliminating eddy-eddy interactions as in the QL approximation (22) of the Boussinesq 
equations. The relative vorticity in the QL simulation is shown in Fig. |B1| The 


positive vorticity anomaly (labeled V) that detaches from a wave crest is initially 
sheared by the mean flow, leading to momentum fluxes that strengthen the mean flow 
and to a decrease of the eddy kinetic energy (Orr mechanism). It is then advected 
around the centre of the cats’ eye (labeled X). But instead of rolling up into a small- 
scale filament as it does in the fully nonlinear simulation (cf. Fig. |^, it moves to the 
western side of the eye, where it joins the wave lobe with positive vorticity west of 
the one from where it originated {T = 5.9). An analogous description applies to the 
negative vorticity anomaly inside the eye. This leads at T = 5.9 to vorticity anomalies 
that have a southeast-northwest tilt of phase lines, consistent with an equatorward 
eddy momentum flux (reflection phase) and an increase of EKE because the mean- 
flow shear goes against the tilt (Orr mechanism). The vorticity map after two cycles 
of absorption and reflection (T = 17.5) is very similar to the initial one {T = 4). 
Differences arise from wave absorption occurring through weakly nonlinear processes 
and, to a lesser extent, from hyperviscosity. This is in sharp contrast to the full 
simulation, in which filamentation takes place, eventually leading to irreversibility (cf. 
r = 4 and T = 17.5 of Fig. [I]). 
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Figure Bl. Evolution of relative vorticity in a QL simulation of the larger- 
amplitude eddies (e = 6) and an Earth-like setting (Ro = 0.06), for which the 
fully nonlinear simulation is shown in Fig.[^ White X’s mark the centres of what 
would become cats’ eyes in the fully nonlinear simulation. The locations labeled 
by V follow a positive vorticity anomaly that detaches from an initial wave lobe. 
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